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The coexistence between different informational molecules has been the preferred mode to cir- 
cumvent the limitation posed by imperfect replication on the amount of information stored by each 
of these molecules. Here we reexamine a classic package model in which distinct information carriers 
or templates are forced to coexist within vesicles, which in turn can proliferate freely through binary 
division. The combined dynamics of vesicles and templates is described by a multitype branching 
process which allows us to write equations for the average number of the different types of vesicles as 
well as for their extinction probabilities. The threshold phenomenon associated with the extinction 
of the vesicle population is studied quantitatively using finite-size scaling techniques. We conclude 
that the resultant coexistence is too frail in the presence of parasites and so confinement of templates 
in vesicles without an explicit mechanism of cooperation does not resolve the information crisis of 
prebiotic evolution. 



PACS numbers; 87.10.-|-e, 87.23.Kg, 02.50.Ey, 02.70.Uu 



I. INTRODUCTION 

The information crisis in prebiotic or chemical evo- 
lution stems essentially from two observations: (i) the 
length of a replicating polymer (i.e., RNA-like template) 
is limited by the replication accuracy per nucleotide |H, 
and (ii) templates that differ significantly from each other 
cannot coexist in a purely competitive setup [2^. Real- 
istic estimates of the error rate of primitive replication 
mechanisms predict a too scanty information content per 
template - less than 100 nucleotides - to permit the com- 
plete codification of the mechanism in just one template. 
Currently an operative replication mechanism requires 
at least three basic sets of different reactions (initiation, 
elongation and termination) Q , thus the primitive infor- 
mation integrator systems must have shared the neces- 
sary information in a number of distinct templates. Yet, 
attainment of template coexistence in a plausible prebi- 
otic scenario is still a highly controversial issue. 

An attractive solution to this crisis is the hypercycle, a 
cyclic reaction scheme in which each replicating polymer 
aids in the replication of the next one, in a regulatory 
cycle closing on itself [4]. This scheme requires that the 
primordial replicators functioned both as templates and 
replicases (i.e., catalysts for replication), a prospect con- 
firmed by the discovery of the catalytic activity of RNA 
in the early 1980s 0, 0]. However, the key assumption 
that each replicator has two separate functions, namely 
a replicase for the next member of the hypercycle and a 
target for the previous member, encountered strong crit- 
icism 0, Q- In fact, natural selection can make each 
element of the hypercycle a better target for replication, 
but it cannot favor the cooperative part of the scheme, 
i.e., to make the replicator a better replicase for other 
replicators. Hence this function is bound to degenerate 
rapidly as natural selection does not protect it against 
deletions and mutations that create the so-called para- 
sites - molecules that do not reciprocate the catalytic 
support they receive. 



Another proposal to resolve the problem of the coex- 
istence between templates, which is in line with the clas- 
sical works on the origin of life f^, [l3| , is to enwrap the 
replicators in isolated compartments or vesicles. It is pre- 
sumed that the information to code for a common repli- 
case is shared among d distinct template types so that 
template replication is feasible only if all template types 
coexist within a vesicle [ll|. In this scheme the repli- 
cators play the template role only whereas the replicase 
role is taken on by a protein - the replicase. In addition, 
it is assumed that a vesicle splits into two daughters af- 
ter a certain number of template copies are produced. 
Alternatively, we may suppose that an effective coupling 
among different template types comes about through a 
common metabolism which is ultimately responsible for 
the survival and reproduction of the vesicle, and that the 
functioning of this metabolism requires the contribution 
of all template types [l2| . 

It should be noted, however, that the proponents of 
the hypercycle have always acknowledged the essential 
function of compartments, particularly in the evaluation 
of the translation products of the information coded in 
the templates 13 1 (see [l3| for the in vitro realization of 
this idea). In reply, the advocates of the so-called pack- 
age models observe that once the replicators are put into 
compartments, the hypercyclic organization is dispens- 
able III. 

So far practically all package models proposed to inves- 
tigate template coexistence (see, e.g., (l^[l5l[l6l[T7l.[T8jl 
assume that the vesicles proliferate or divide with a rate 
that depends on their template compositions. This as- 
sumption can be interpreted as a group selection pressure 
acting at the vesicle level to favor vesicles of a particu- 
lar makeup. Despite years of intensive research on vesi- 
cle dynamics, however, there is no experimental evidence 
that the vesicle fission rate could depend on the nature 
of the chemicals confined inside it jlj,, 20] . Hence a more 
conservative stand is to admit that vesicle fission is trig- 
gered by the total concentration of the confined templates 
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rather than by their individual proportions. (Of course, 
the total template concentration does depend on the ex- 
istence of the replicase and hence on the presence of all 
d template types.) Interestingly, in their seminal work 
Niesert et al. take this cautious position, but to avoid 
the unbounded growth of the vesicle population, they 
discard supernumerary vesicles according to an arbitrary 
prospective value which essentially gauges the odds of a 
vesicle to leave viable descendents Then the result- 
ing model becomes very similar to the group selection 
models mentioned above. 

In this paper we propose and study analytically a vari- 
ant of the original package model of Niesert et al. where 
the number of vesicles is unbounded and no fitness or 
prospective values are assigned to the vesicles. In par- 
ticular, we derive a recursion equation for the average 
number of vesicles with a given template composition, 
and a set of equations for the extinction probability of 
the distinct vesicle types. Our results indicate that an 
important conclusion of the original work - high values 
of the replicase processivity compromises the viability of 
the population in the presence of parasites - is proba- 
bly an artifact of sampling the vesicles according to a 
prospective value. In this line, we use finite size scal- 
ing to show how a strategy of discarding supernumerary 
vesicles at random can efficiently recover the analytical 
results. 



II. MODEL 

We follow the original package model proposed by 
Niesert et al. [Tlj and consider a metapopulation com- 
posed of a variable number of vesicles, each of which en- 
closes a certain number of templates. There are d dis- 
tinct functional types of templates / = 1, . . . , d and a 
non-functional type I — 0, termed parasite, which has 
an impaired function but an unchanged replication rate. 
Due to imperfect replication, functional templates mu- 
tate to parasites with probability u. Back mutations as 
well as mutations between functional templates are ne- 
glected. To be consistent with the conjecture that all 
templates display identical targets to the replicase since 
they derive from a common ancestor (most likely were 
members of a same quasispecies [ij), the replication rate 
is assumed to be the same for all templates (including 
the parasite). Here we do not contemplate two additional 
processes allowed for in the original model, namely, the 
possibility of mutation to lethal genes or the possibility 
of accidents. Both actions prompt the immediate demise 
of the vesicle. 

The life cycle (i.e., one generation) of the metapop- 
ulation comprises three events - template replication, 
vesicle fission and vesicle extinction - that take place in 
this order. In this contribution we modify the first two 
events in order to produce an analytical formulation of 
the metapopulation dynamics. In particular, we assume 
that the quantity of templates confined in each vesicle 



before replication is fixed to a certain value A. Tem- 
plate replication doubles this number but then fission of 
the mother vesicle into two daughters of identical size re- 
stores it to the original value. Hence A can be interpreted 
as the number of replicated molecules between two vesicle 
fissions, which is essentially the processivity of the repli- 
case, i.e., the number of template copies the replicase can 
produce in a unit of time, taken here as the time between 
two consecutive fissions. 

In contrast, in the formulation of Niesert et al. [ll| 
the two daughter vesicles can have different sizes s = 
0, . . . , S" and S — s, with s distributed by the binomial 

distribution | | 2^"^ where S is the size of the mother 

vesicle after template replication, so that the number of 
templates can vary among the vesicles. The processivity 
A of the replicase, however, is the same for all vesicles 
and so, in the average, the number of templates within 
each vesicle equals A. Thus, essentially, our formulation 
neglects fluctuations in the number of templates inside 
the vesicles. As we will show in Sec. |Vl these variations 
in the modeling of the vesicle dynamics do not change 
qualitatively the main results of the model. 

In our model, the composition of each vesicle is fully 
characterized by the vector k = [ko, . . . , kd), where 
the entries fc; yield the number of templates of type 
I = 0, . . . , d in the vesicle and satisfy the constraint 

J2f=o^i — ^- So there are exactly A/r — ^^^"^^ 

distinct types of vesicles - the number of compositions 
of A into d -I- 1 parts. To keep track of all vesicle types 
we use a combinatorial algorithm to generate and label 
those compositions plj . 

The more restrictive assumption of the model is prob- 
ably the choice of equal replication rates for the distinct 
types of functional templates as well as for the parasites. 
The supposition that the parasite and the functional 
classes have equal replication rates is plausible since a 
parasite is essentially a functional template whose activ- 
ity was impaired by a mutation in the region coding for a 
piece of the replicase. In any event, allowing the parasites 
to replicate faster or slower than the functional templates 
has an effect similar to that of increasing or decreasing 
the mutation probability u. The choice of different repli- 
cation rates for the functional templates, however, has 
drastic consequences in the limit of large A, where the 
intra-vesicle dynamics becomes deterministic: the more 
efficient template type drives the other functional tem- 
plates to extinction, thus preventing coexistence even in 
the absence of parasites. In fact, in a class of models 
where the number of vesicles is fixed and very large, there 
is a limiting value of A above which template coexistence 
is impossible [f^, [3- We note, however, that a more 
realistic scenario would allow the replication rates of the 
functional templates to vary under the pressure of nat- 
ural selection. Since only the exact balancing of those 
rates guarantees coexistence (and so survival) for large 
A, one expects the selection of this ideal symmetric set- 
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ting. This reasoning supports the assumption of equal 
replication rates for the functional templates. 



A. Template replication 

We assume that the replication of the A templates 
encaged in a vesicle follows a Wright-Fisher process in 
which the A offspring are chosen in parallel ^] . Admit- 
ting equal replication rates and unidirectional mutation 
to the parasite class, the probability that a set of tem- 
plates ko, . . . ,kd produces the set of offspring io, . . . , is 
given by the multinomial distribution 

i? (? I fc) = . . . [wo + uil- wo)Y" n l^i (1 - ■ 

(1) 

where wi = A://A for / = 0, . . . , d so that J2i=q^i — 1- 
The interpretation of Eq. ([T]) is straightforward. On the 
one hand, the probability of producing a functional off- 
spring of type I is given by the probability of choosing 
a template of the same type, wi, times the probabil- 
ity that the copy produced is faithful, 1 — u. Parasites, 
on the other hand, are produced by unfaithful copies of 
functional templates with probability u {wi + . . . + w^) 
or by copies of parasites themselves, with probability 
Wo- Once the template replication process is completed, 
we are left with a vesicle of size 2A and composition 
ka + io, . . . ,kd + id- This procedure is repeated for all 
vesicles in the metapopulation. 

In the original model [ll| the replication procedure is 
sequential rather than parallel. For a given vesicle we 
choose a template at random (from those inside the vesi- 
cle) and make a copy of it. If the template is of a func- 
tional type then the copy will become a parasite with 
probability u. If the chosen template is a parasite then 
the copy will also be a parasite. Both template and copy 
(corrupted or not) are returned to the vesicle and the 
process is repeated A times, so exactly A new templates 
are added to the vesicle. This difference in the modeling 
of the template replication process does not affect the 
results in any significant way. 



B. Vesicle fission 

The doubling of the size of the vesicles caused by the 
template replication process leads to the splitting of the 
vesicle in two daughters. Our simplifying assumption 
here is that the vesicle of size 2A splits into two vesicles 
of size A. The assignment of the A templates to one of 
the daughter vesicles is modeled by a process of sampling 
without replacement which is described by a multivariate 
hypergeometric distribution. Explicitly, given the com- 
position of the mother vesicle after template replication 
k + i, the probability that one of the daughter vesicles 



has composition mo, . . . , rud is simply 




with mi < ki + ii and X]i=o ™' ~ ^- course, if one of 
the daughter vesicles is described by m, then the other 
will be described by fc + i — m. The random assortment of 
templates to the daughter vesicles is the only mechanism 
responsible for the loss of the essential genes for survivor- 
ship, a phenomenon termed assortment load. The loss of 
a functional template occurs when it is assigned to an 
inviable vesicle, i.e., a vesicle that does not contain the 
complete set of functional templates. 

As mentioned before, in the original model [ll| the 
sizes of the daughter vesicles are binomially distributed 
random variables and so some vesicles can become very 
large since what prompts vesicle fission is not its absolute 
size, but the production of A template offspring. This 
asymmetry in the fission process renders the population 
more susceptible to the presence of parasites (sec Sec. 
W} . but as already said, does not change qualitatively 
the results of the model. 



C. Vesicle extinction 

The viability of a daughter vesicle is guaranteed pro- 
vided it has at least one copy of each functional tem- 
plate. Any vesicle lacking one of those templates is dis- 
missed. Strictly, we do not need to assume that the in- 
viable vesicles disappear from the metapopulation, but 
since the templates caged in those vesicles are unable 
to replicate - their replicase is not codified for - there 
is no point to follow their evolution any further. We 
note that the total number of viable vesicles A/y = 
X]feo>o fei>i ^ (X]f where 6{m, n) is the Kronecker 



delta, is simply 




. To the leading order in A both 



quantities JVy and JVt increase as A'^ and for large A we 
find Afv/AfT = l-(f/k + (A-2) . 



D. Metapopulation dynamics 

^From the processes described above, it is clear that 
the size of the metapopulation (i.e., the number of viable 
vesicles) can, in some cases, increase without bounds. 
Such unbounded growth renders a direct simulation ap- 
proach of the vesicle population dynamics unfeasible, ex- 
cept for the few initial generations. To circumvent this 
difficulty, here we derive a set of recursion equations for 
the average number of vesicles of type rfi at generation 
denoted by $t (m). 
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The basic idea is to derive a transition matrix that con- 
nects the mother vesicle k with the two daughter vesicles 
rfia and m;,. From Eq. ([2|) we can immediately write 
down the transition probability from the mother vesicle 
to the first daughter, 



Ga (ma I fc) =^^^(rn, | fc + ?) i? | k 



(3) 



The derivation of the transition probability from k to 
the second daughter rhb = k + i — rha is more involved 
because of the dependence on the intermediate states i 
which we ultimately want to sum over, as done in Eq. 

Given that the first daughter has composition rua, 
the probability that the second daughter has composition 
rfib is 



H {rhb I ma) 



F [rha I mb + rha) 



G 



a {rha I fc) 



R ( rhb 



k I k 



J4) 

which is obtained by considering only the term i = rhb + 
rha — k in Eq. ([3]) and properly normalizing. Clearly, the 
joint probability that the daughter vesicles are of types 
rha and rhb given that the mother vesicle is of type k is 
simply 

Pab (rha, rhb\k^ ^ H {rhb I rha) Ga (rha I k^ , (5) 
and so the desired transition probability is 
Gb (rhb I ^) = X! (rha, rhb \ kj 

= F (^i I m^, + i? -\- i — k \ ^ 

i 

(6) 

where we have replaced the dummy index rha by i to 
facilitate the comparison with Eq. The transition 

matrices given by Eqs. ([3]) and ([6]) allow us to write a 
recursion equation for the average number of the different 
types of vesicles in the metapopulation. 



t+i 



{rh)=J2'^t (k 



Ga 



k]+Gi 



(rh I kj 



(7) 

where the primed sum is over viable vesicles, i.e., vesicles 
that contain at least one copy of each functional gene, 
ki > for / > 0. This restriction is the expression of the 
vesicle extinction process. 

In principle, the solution of the recursion equations 
([7]) yields detailed information on the time evolution of 
the metapopulation. But the computational resources 
needed to generate the entries of the matrices Ga and 
Gb seriously constrain the range of A and d that can be 
studied in practice. In this contribution we show how 
this difficulty can be circumvented by considering a fi- 
nite population of vesicles with the same growth rate 



per generation as the ideal, unrestricted metapopulation 
described above. Nevertheless, some interesting infor- 
mation can be obtained from the analytical approach as 
described in Sec. IIIII 



E. Extinction probability 

To point up the stochastic nature of the underlying 
vesicle dynamics - the approach based on the recursion 
equations ([7]) is deterministic as the focus is on the aver- 
age number of a given vesicle type - here we describe a 
general formulation to calculate the extinction probabil- 
ity Pe (itj of the lineage sprouted by a vesicle of type k. 

Generalizing the classic approach to evaluate the extinc- 
tion probability in the Galton- Watson process [H, we 
can write the following set of equations 



Pr k 



Pab (rha,rhb \ k) Pe {rha) Pe {rhb) (8) 



with the convention that Pg {rh) = 1 if the vesicle of type 
rh is inviable. Note that P^ (k^ = 1,VA: is always a solu- 
tion. Surprisingly, this system of A/V nonlinear coupled 
equations easily yields to the simple iterative solution 

method that begins with the guess Pe (itj = for all 
viable vesicles. 

III. UNRESTRICTED GROWTH 

It is clear from Eq. ([7]) that the asymptotic regime of 
the dynamics is characterized either by the simultaneous 
divergence or by the simultaneous vanishing of $oo {rh) 
for all viable vesicles. The main goal here is to find the 
values of the control parameters A, d and u that separates 
these two regimes. In general, this critical parameter 
setting can be found by direct numerical iteration of the 
recursion equations. 

Let us discuss first the simpler, extreme case A = d 
for which there is only one type of viable vesicle, namely, 
m* — (0, 1, 1, ... , 1). In this case, Eq. ^ simphfies con- 
siderably and allows us to carry out analytically the sum- 
mations in Eqs. ^ and We find Ga {rh^ \ m*) = 
Gb {rh^ I rh^) = Qd{u) where 



ild {u) 



d 

E 

i=0 



(1-^)7 



2d 
d 



(9) 



and so ^t+i{rh^,) = 2fld {u) {rh^,). A straightfor- 
ward numerical evaluation of Eq. ^ for u — yields 
^2 (0) — 7/12 and il.d>2 (0) < 1/2. Since is a mono- 
tonically decreasing function of u, the latter inequality 
implies that ^d>2 {u) < 1/2 for nonzero u as well. This 
indicates that template coexistence is unattainable for 
A = d > 2. For A = d = 2, however, the picture is 
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different: the average size of the vesicle population will 
increase exponentially with increasing t provided u < Uc 
where Uc = 3 — 2\/2 is the solution of ^2 {u) — 1/2. We 
find this simple analytical result reassuring because it 
proves that, even for finite vesicle capacities, functional 
templates can persist in the presence of a steady drain 
towards the parasite class. 

The evaluation of the extinction probability is also 
straightforward in the case A = d. Using Eq. (O we 
write the probability that a viable vesicle produces two 
viable daughters (i.e., vesicles of type m*) as 



P2 = d\ 



2^ 



2d 

d 



(10) 



Now, Eq. ^ yields the probability that the first daugh- 
ter of a viable vesicle is also viable, regardless of the con- 
dition of the second daughter, so the probability that 
a viable vesicle produces a single viable offspring, no 
matter whether it is the first or the second daughter, 
is pi = 2{Qd —P2)- Hence the probability of produc- 
ing two inviable daughters is po = 1 — pi — p2- In this 
case Eq. reduces to the simple quadratic equation 
Pe = Pa +PiPe +P2Pe with = Pe ("I,). The two so- 
lutions are Pe = 1 and Pe = {2Ud — 1) /p2. The latter is 
physical provided > 1/2 which, as pointed out before, 
holds only for d — 2 and u < 3 — 2V2. 

For A > d, we have to resort to the numerical iteration 
of Eq. (O or to the numerical solution of Eq. ([5]) to obtain 
the critical parameter setting that determines the regime 
of viability of the metapopulation. In the former method, 
we begin the iteration (t — 0) with a single parasite-free 
vesicle, whose composition of functional templates is as 
balanced as possible, e.g., (0, A/d, . . . , A/d) in the case of 
integer A/d. 

The results for the case of perfect replication accuracy 
M = are summarized in Fig. [TJ The critical processivity 
value Ac above which the size of the metapopulation di- 
verges is very well described by the fitting Ac = d^/2 as 
shown in the figure. This indicates that the assortment 
load can be compensated for if the redundancy A/d is 
larger than half the diversity value, i.e., provided that 
each vesicle contains at least d/2 copies of each func- 
tional template. This simple result shows that there is 
no fundamental impediment to the coexistence of an ar- 
bitrary number of template types in the case of error-free 
replication if the cost of redundancy is neglected. To un- 
derstand the scaling Ac ^ d^ at the critical boundary in 
the error-free replication limit we must look at the ratio 

A-1 



r between the number of viable vesicles 



d- 1 



and 



the total number of vesicles 



A 



-d- 1 
A 



, given by 



i/A 



+ z/A' 



(11) 



We note that since the parasite class is not taken into ac- 
count in this error-free replication analysis, the number 
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FIG. 1: Logarithmic plot of the value of the processivity Ac 
above which the population size diverges against the number 
of functional templates d in the case of error-free replication 
u = 0. Below Ac the extinction of the lineage is certain. The 
solid line is the fitting Ac — 
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FIG. 2: Mutation probability Uc above which the metapop- 
ulation is inviable as a function of the replicase processivity 
A for template diversity d as indicated in the figure. The 
lines are the analytical results for unrestrained growth and 
the symbols are the results of the finite-size scaling analysis. 



of viable vesicles as well as the total number of vesicles 
differ from the quantities My and Nt introduced before. 
The only way to obtain nontrivial values of this ratio (i.e., 
r ^ 0, 1) for large A and d is to suppose that d^/A re- 
mains of order of 1. In this case we find r ~ exp (— d^/A) 
and so = e^^ « 0.135 at the critical boundary. 

We turn now to the analysis of the case where par- 
asites are allowed, i.e., m > 0. Fig. [2] summarizes the 
main results, namely, the dependence on A and d of the 
critical mutation probability Ue above which the lineage 
is inviable. The curves intersect the axis Uc = at the 
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FIG. 3: Survival probability of the lineage produced by a sin- 
gle balanced, parasite-free vesicle as function of the mutation 
probability u for d = 2 and (left to right) A = 2, 3, 4, 5, 6, 8, 10, 
and 20. The lines are the analytical results for unrestrained 
growth and the symbols are the results of the finite-population 
simulations with A'^ = 10^ and 10^ independent samples. 

values of A exhibited in Fig. [1] The remarkable result 
revealed in Fig. [5] is that, for fixed d, there exists a value 
of the mutation probability above which coexistence be- 
tween the d functional templates is impossible regardless 
of the replicase processivity value or, equivalently, the 
redundancy value. This result follows from the fact that 
Uc tends to a well-defined value less than 1 in the limit 
A — !■ oo. This is reminiscent of the error threshold tran- 
sition of the quasispecies model for which the replication 
fidelity limits the length of the templates and hence the 
amount of information that can be stored in the molecu- 
lar population [Ij] . Here the limitation is on the number 
of different types of functional templates that can coexist 
within a vesicle and so on the total amount of information 
that can be stored in the vesicle. Another result shown in 
Fig. [5] is the impracticability of the analytical approach 
for large N't- the time required to evaluate the matrix 
entries defined in Eqs. ([3]) and ^ is simply prohibitive so 
the curves are truncated at the values of A that surpass 
our computational resources. Fortunately, the analysis 
of a finite population can greatly extend these limits, as 
we will show in Sec. IIVI 

In addition to the threshold values exhibited in the pre- 
vious figures, the analytical approach allows us to obtain 
some detailed information about the composition of the 
metapopulation and the nature of the stochastic process. 
In particular, in Fig. [3] we show the survival probabil- 
ity Ps = 1 — Pe oi the lineage sprouted by a balanced, 
parasite-free vesicle obtained by solving numerically Eq. 
([8]) for d = 2. For A = 3 and 5 the ancestor vesicle is 
(O, ^^^)- Since the template dynamics becomes de- 
terministic in the limit A — > oo, the survival probability 
must tend to a step function as indicated in the figure. 

Interestingly, we find that regardless of whether the 



FIG. 4: Fraction of functional templates per viable vesicle 
in the steady-state regime for d = 2 and (top to bottom) 
A = 3, 4, 5, 6, 8, 10, 20, and 40. The symbol x indicates the 
value of u beyond which the population is inviable. 

process is subcritical (u > Uc) or supercritical (u < Uc) 
the fraction / of functional templates per viable vesicle, 
and consequently the fraction of parasites, rapidly tends 
to a steady-state value. This fraction, defined by 

f = hm ^ ^ (12) 

AE£$*(fc) 

is shown in Fig.|3]as a function of the mutation probabil- 
ity. For large A we find that / vanishes as A~^ provided 
u is nonzero and, in particular. A/ = d = 2 for u = I. 
So in this limit there is only a finite number of functional 
templates in each vesicle, whereas the number of para- 
sites grows linearly with A. The reason the population 
is viable (see Fig. [3]) even in these circumstances is that 
the chances for choosing functional templates for repli- 
cation in some vesicle is not negligible when the number 
of vesicles is greater than A, which is always the case in 
the supercritical regime after a few generations. 

IV. FINITE POPULATION 

As pointed out before, the impracticability of generat- 
ing the entries of the matrices that govern the transitions 
between viable vesicles for large values of A and d limits 
the applicability of the analytical solution summed up in 
Eqs. (O and ([8]). To get around this obstacle we consider 
here an alternative approach based on the Monte Carlo 
simulation of a finite population. In contrast to classical 
models of populations genetics (e.g., the Wright-Fisher 
and Moran models (2^ 1 in which the population size is 
kept fixed, here we allow the number of vesicles to vary 
from to a fixed maximum value N. The idea is to 
implement the dynamics exactly as done in the case of 
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FIG. 5: The survival probability of a population composed 
initially by m = balanced, parasite-free vesicles for K = d = 
2 and (left to right) iV = 40, 80, 100, 200 and 1000. The solid 
line is the analytical result for unrestrained growth. Each 
symbol represents an average over 10^ samples. 

unrestrained growth, except that whenever the number 
of vesicles becomes greater than iV, the surplus vesicles 
are discarded randomly. The discard takes place before 
the check of the viability of the vesicles (extinction pro- 
cess). This scheme is reminiscent of the so-called Russian 
Roulette used in the Monte Carlo simulation of neutron 
production in nuclear reactors [25j . 

In the subcritical regime, the introduction of the up- 
per bound N is innocuous since the population size is 
likely to remain small before the extinction outcome any- 
way. In the supercritical regime, however, a too small 
bound may prevent the lineage to produce and retain 
a minimum number of viable vesicles that would avoid 
extinction and so one expects Uc{N) < Uc- (An opera- 
tional definition of Uc{N) will be given later.) The finite 
population scheme is effective in the practical situation 
N ^ A/y provided that only a small fraction of the Nv 
viable vesicles would actually be present in the metapop- 
ulation if it were allowed to grow unrestrained. 

The previous setup for the initial structure of the pop- 
ulation - a single balanced, parasite-free vesicle - is not 
suited to study finite size effects on the estimate of the 
critical mutation probability, because there is no oper- 
ational way to define Uc{N). For that end an effective 
strategy is to begin with a population of m identical such 
vesicles. Since the vesicles evolve independently there is 
a simple relationship between the probability that a pop- 
ulation with initial size m thrives, denoted by PJ", and 
the survival probability of a single vesicle Ps exhibited in 
Fig. [31 namely, 

P;" = l-(1-P,)'". (13) 

For m — > cx) this quantity tends to a step function that 
takes on the values 1 if u < Mc and otherwise. In the 



FIG. 6: Mutation probability Uc{N) at which 50% of samples 
of a population of A'^ balanced, parasite-free vesicles survive 
for d — 2 and A — 2, 3, and 4. The linear fittings (dashed 
lines) Uc{N) = Uc — ua/N allow us to estimate Uc- 



finite population simulations we set m — N, since our 
focus in on the behavior of P™ when both quantities - 
the initial size m and the size upper limit N - become 
arbitrarily large. The results for Pj^ are shown in Fig. 
[5] for different population sizes N. As N increases, the 
finite-population results approximate those for the un- 
restrained growth represented by the step function. To 
quantify this approach, we arbitrarily define Uc (N) as 
the value of the mutation probability at which = 1/2 
so that the critical value Uc for — s- oo can be inferred 
as illustrated in Fig. [HI 

Use of P^ instead of Pg is crucial for this analysis, 
since regardless of the definition of Uc{N) (e.g., we could 
define it as the mutation probability at which P^ = x 
for any < a; < 1) this quantity tends to Uc in the limit 
of A^ large. The extrapolated values to 1/A^ — > are 
presented in Fig. [2] and agree perfectly with the available 
analytical predictions. This gives us confidence to use 
the finite population estimates in the cases where the 
analytical approach is not practical. 

Fig. [7] summarizes the results of the data collapsing 
method (see, e.g., [H, |23|) applied to the data of Fig. 
[5l The survival probabilities Pj^ collapse into a sin- 
gle universal form (scaling function) if plotted against 
the scaled mutation probability {u — Uc) N where Uc is 
the critical mutation probability of the infinite popula- 
tion. This scaling function shows that, for u — Uc fixed, 
liniAT^oo Pj^ tends to 1 if u < Uc and to otherwise, and 
that the characteristics of the threshold transition per- 
sist across a range of u of order 1/N about Uc- A similar 
finite-size scaling analysis was employed to fully charac- 
terize the error threshold transition of the quasispecies 
model in the cases where the population size as well as 
the molecules lengths are fixed and finite [28j . 

Of course, the finite population scheme can be used to 
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FIG. 7: The data of Fig. [5] (d = A = 2) plotted against the 
scaled mutation probability {u — Uc) N with Uc = 3 — 2\/2. 
The collapse of the data into a single A'^-independent function 
signals the occurrence of a threshold phenomenon at Uc- 



0.3 




FIG. 8: Critical mutation probability obtained via finite-size 
scaling for the package model proposed by Niesert et al. (A) 
and the variant proposed in this paper (Q) for d = 2 and 3 
as indicated. 



study the original setup in which the initial population 
comprises a single vesicle, m = 1, as well. In this case the 
size limit N = 1000 suffices to obtain perfect agreement 
with the analytical results, as shown in Fig. [31 



V. DISCUSSION 

In Sec. |TT]we have outlined the differences between the 
original package model proposed by Niesert et al. and 
our more tractable variant. Here we present a brief com- 
parison between the main predictions of these models. 
The first important observation is that in the absence of 



parasites (w = 0) both models yield identical results for 
the viability boundary in the plane (d. A) (see Fig. [T]). 
When parasites are present, however, there is a substan- 
tial quantitative difference between the critical mutation 
probability of the two models, as shown in Fig. [S] The 
scheme based on synchronous template replication and 
symmetric fission seems to be considerably more robust 
to the action of parasites than the less structured proce- 
dures of the original proposal. This result supports the 
view that the mechanisms of segregation of modern cells 
originated in response to mutation pressure [29| . 

It is hard to see why parasites are more harmful in the 
asynchronous replication and asymmetric fission setting. 
Considering that parasites are rare at the beginning, their 
spread should be hampered by the asynchronous replica- 
tion scheme in the initial generations and then speeded 
up when the parasites become numerous. By simulating 
the two template replication schemes with the same fis- 
sion mechanism, we have verified that the choice of the 
form of update - parallel or sequential - practically does 
not affect the critical mutation probability Uc- Thus the 
key element to explain the quantitative differences be- 
tween the models illustrated in Fig. [5] must be the fission 
mechanism. To get some insight on that, let us consider 
the situation where a mother vesicle of size 2A contains 
two functional templates of a certain type. Clearly, from 
the mother vesicle's perspective the optimal strategy is 
to send one template to each of her daughters. The prob- 
ability this happens for the asymmetric fission strategy is 
1/2 independently of the vesicle size. The symmetric fis- 
sion scheme in turn yields a slightly larger probability for 
this event, namely, 1/2 x (1— 1/2A)~ . This tendency 
of the symmetric fission strategy to a more balanced dis- 
tribution of templates of the same type to the daughters 
is probably the reason of its enhanced robustness against 
parasites. 

Our finding that Uc is a nondecreasing function of A 
(see Fig. [5]) is at variance with the results of Niesert et 
al. which predict that Uc would reach a maximum and 
then decrease towards zero as A increases further llj. 
The reason may be the criterion for discard of supernu- 
merary vesicles used in that work, which was based on 
three properties: the degree of equipartition of the copies 
among the different functional templates, the number of 
parasites and the overall redundancy of the functional 
templates. In fact, we have verified (see Fig. S]) that for 
large A and not too low u, the surviving vesicles in the 
supercritical regime are heavily loaded with parasites and 
so use of such selection criterion would purge them from 
the population resulting in a premature extinction. 

Although the finite population simulations were used 
here as a tool to validate and complement the analytical 
results, they are of interest on their own. In particular, 
the Muller ratchet [13, HH and the mutational meltdown 
[3^ are important stochastic phenomena that result in 
the accumulation of mutations in finite populations (see 
[s^ l for the study of both phenomena in growing lin- 
eages). In our framework, the counterpart of accumu- 
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lation of mutations is the accumulation of inviable vesi- 
cles, which is explicitly ruled out by the assumption that 
those vesicles are unable to divide into daughter vesicles. 
This peculiar aspect of the model was severely criticized 
by Eigen et al. [13] who pointed out that the vesicles in 
the model Niesert et al. [ll| cannot evolve because of 
that assumption. Since there is no competition among 
the vesicles in the case of unrestricted growth, allowing 
the inviable vesicles to divide as well would have no effect 
at all on the dynamics of the viable vesicles because it is 
not possible to produce a viable vesicle by fissioning an 
inviable one. In the finite population case, on the other 
hand, the inviable vesicles would accumulate steadily and 
ultimately would reach fixation in the metapopulation. 

VI. CONCLUSION 

The goal of the research on prebiotic evolution is to 
put forth a coherent scenario for the origin and early de- 
velopment of life. So at this stage it is appropriate to 
appraise the main results of our analysis of this classic 
package model, summarized in Fig. [21 Given the sponta- 
neous error rate per nucleotide e and the molecule length 
L ^ 1 we can readily obtain the value of the probabil- 
ity of mutation from functional templates to parasites, 
u = I — exp(— eL). A plausible estimate for these pri- 
mary parameters is e ~ 10^^ and L ~ 100 [1] which 
yields u ^ 0.6. A glance at Fig. [2] leads to the disas- 
trous conclusion that even the coexistence between two 



templates is prohibited in these circumstances. It is in- 
structive also to compare our results with those of the 
hypercycle which guarantees the stable coexistence of at 
most d = 4 templates Q (see [sJl for the analysis of the 
hypercycle in the presence of an error tail class similar 
to the parasite class considered here). According to Fig. 
[21 rf = 4 functional templates can coexist provided that 
u < 0.3 which implies that L < 35, resulting in a total 
of 140 nucleotides, a meager improvement over the 100 
nucleotides prediction for a free replicator. 

It should be observed, in addition, that the critical 
mutation probabilities Uc exhibited in Figs. [21 and [Hare 
best case results since neither lethal mutations nor acci- 
dents were considered in our calculations. Hence, con- 
trary to the claims of Niesert et al. [ll[, the kind of 
template coexistence achieved in a simple package model 
does not resolve the prebiotic information crisis. Special 
mechanisms to prevent independent information carriers 
from competing with one another within the compart- 
ment must be posited. Although unwarranted assump- 
tions like the hypercyclic organization Q or the cou- 
pling between the templates and the package metabolism 
may weaken the credibility of the models, 
so far they seem to stand as the only options to tackle 
the information crisis of prebiotic evolution. 
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